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1 The Problem Considered 

The present numerical study is concerned with the fundamental physics of the multi-way interaction between 
turbulence, chemical reaction and buoyancy in a nonpremixed flame. The method of direct numerical 
simulation (DNS) is used to solve the instantaneous three-dimensional governing equations. Because of the 
present supercomputer limitations, we consider two simple flow geometries, namely an initially uniform flow 
without shear ( equivalent to grid-generated turbulence) and an initially uniform shear flow. In each flow, 
the fuel and oxidant initially exist as two separate streams. As the reactants mix, chemical reaction takes 
place and exothermic energy is released causing variations in density. In the presence of a gravity field, the 
spatial and temporal distributions of the induced buoyancy forces depend on the local density gradients and 
the direction of the gravitational acceleration. The effects of buoyancy include the generation of local shear, 
baroclinic production or destruction of vorticity, and countergradient heat and mass transport. Increased 
vorticity and small-scale turbulence promote further mixing and reaction. However, if the strain-rates 
become too high, local flame extinction can occur. 

Our objective is to gain an understanding of the complex interactions between the physical phenomena 
involved, with particular attention to the effects of buoyancy on the turbulence structure, flame behavior 
and factors influencing flame extinction. 


2 Mathematical Description 
2.1 Approach 

We focus our attention to a location in the vicinity of the interface between two initially homogeneous 
gaseous reactant streams in an unbounded domain. The three-dimensional flow field is initially homogeneous 
isotropic or sheared turbulence with subsequent development due to density variation arising from chemical 
energy release and gravitational effects. A schematic of the flow configurations considered is shown in 
Figure 1. A simple one-step irreversible chemical reaction is considered: 

v f F + voO — ► i/pP , ( 1 ) 

in which two nonpremixed reactants (fuel F and oxidant O) react to form product P. Both finite and infinite 
reaction rates are considered. In the former, the reaction zone has finite thickness (flamelet), whereas 
the latter corresponds to the flame sheet approximation. The finite reaction rate is a function of both 
concentration and temperature. In the case of an infinite reaction rate, a conserved scalar formulation is 
employed. Transport properties of the fluid such as the mixture dynamic viscosity //, thermal conductivity 
I \ , and mass diffusion coefficient pV are assumed constant for the initial stages of the work. The dependence 
of these properties on temperature and concentration will be considered later in the study. 
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2.2 Governing Equations 

We consider a low Mach-number turbulent flow with density variations arising from chemical energy release. 
Thus, the kinetic energy is small in comparison to the thermal energy. Since the time scales of the high 
frequency acoustic waves are several orders of magnitude smaller than those associated with the convection 
processes, the former do not interact effectively with the flow dynamics, and thus their effects can be filtered 
out from the fully compressible equations [6]. The resulting governing equations in nondimensional form are: 


^ + V-(/m)= 0, 

(2) 
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9 + V • (puu) = - Vpd) - (1 /Re)V ■ r + ( pe k )/Fr , 

at 

(4) 

p 9 !. + p(„ . V)T = (1/Pr Re) V • (VT) + Da Ce r T , 

at 

(5) 

9 (pY F ) + v . ^Yp) - (1 /Sc Re) V • (VYf) + Da r F , 

at 

(6) 

d(pYo) + v (/w y o) _ (1/5c Re) v . (V y o) + D ar 0i 

at 

(7) 

P= P ( 0 ) /T, 

(8) 

where standard nomenclature is used and e * is a unit vector in the gravity direction. 

The pressures 


and pi 1 ) are the (zeroth order) hydrostatic (or thermodynamic) pressure and the (first order) hydrodynamic 
pressure associated with fluid motion, respectively. In an open domain as considered here p (0) is constant. 
We also split p f 1 ) into two parts, one is independent of gravity effects and the other equals (p 0 g z), thus we 
have in dimensionless form: Vpl 1 ) = VP + e k / Fr. 

Equation (4) then becomes: 

^ (pu) + V ■ (puu) = -VP-(l/Re)V -T + [(p-l)e k )\/Fr. (9) 

The non-dimensional reaction rates are 

r T = ( pY F y F (pY 0 y° e~ Ze ' T , (10) 


, i= m^ rr 

Vvp up 


(ii) 


where Wi is the molecular weight of species i and v[, v" are the stoichiometric coefficients for species i existing 
as reactant and product, respectively. The Reynolds number, Froude number, Prandtl, and Schmidt numbers 
are defined as 


Re = PoUoLo/fio , Fr = U 0 2 /(L 0 g) , 


Pr = fi 0 C Po /K 0 , Sc = /i 0 /(p 0 ’P 0 ) 5 (12) 

The subscript ’o’ denotes the reference value for the particular quantity. The Damkohler number, defined 
as the ratio of the convective time to a characteristic chemical reaction time, is given by: 

Da = (L 0 /U o )/ B , (13) 

where B is the pre-exponential coefficient in the dimensional expression for tt : 

r T = (pY F /W F y F ( pY 0 /WoY 0 B e ~ E ° ' . (14) 

The non-dimensional energy release parameter and activation energy (Zeldovich number) are: 

Ce = [hpWF^F + hoWovo — h < pWpvp]l\WFVFC Po T 0 } , 

Ze = E 0 /RT 0 . (15) 
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In the limit of an infinite reaction rate, the reacting mixture is effectively at chemical equilibrium and the 
temperature and composition can be determined from a single conserved scalar. In this case, the above 
energy and species concentration equations are replaced with a single equation for the mixture fraction [8]. 

2,3 Numerical Procedure 

The computational domain considered in the simulations is a finite cube of side length =1. The governing 
equations are discretized in space using finite differencing except those terms corresponding to the mean 
advection where a pseudospectral (Fourier) interpolation is used. In general, a finite-difference method 
requires a larger number of grid points for a given resolution as compared with a fully spectral method. 
However, the former allows more flexibility in prescribing the boundary conditions. Spatial homogeneity in 
the horizontal plane allows periodic boundary conditions to be imposed in the x— and y— directions. In 
the vertical direction, z, a convective outflow boundary condition is employed in order to allow expansion 
of the flow due to heat release. The simulations are initialized with divergence free velocity fields having 
random (Gaussian) fluctuations of prescribed correlation spectra. The initial scalar distribution represents 
two separate streams of (nonpremixed) reactants. 

The basic numerical integration scheme follows that of [4]. Integration in time is performed via the Adams- 
Bashforth scheme. Pressure is treated implicitly and obtained using a Poisson solver. The basic procedure 
has been used to simulate both decaying isotropic and homogeneous sheared incompressible turbulence ([1] 
[2] [4] [3] [7]), and an exothermic turbulent nonpremixed reacting flow with zero gravity and infinite rate 
chemistry [8]. In the present work, further modifications have been made to include finite rate chemistry 
and nonzero gravity. 


3 Results 

Due to the imposed 6-page limit, we present only the results for Da = oo but with different values of Fr 
to examine the effects of gravity. As discussed earlier, for the case of Da = oo we replace equations (5), 
(6) and (7) with a single equation for the conserved mixture fraction. Distributions of the temperature and 
mass fractions of all the species are then obtained algebraically from that of the mixture fraction [8] . The 
results are obtained from a 96 3 mesh and initially homogeneous isotropic turbulence with Re\ y0 = 25, and 
the nondimensional energy release parameter, Ce, is prescribed such that the ratio Tfi am e/T 0 = 6. These 
conditions insure the resolution of the Kolmogorov scale motion throughout the simulation. No mean shear 
is imposed on the flow. First we discuss the main features of the flow in zero gravity (Fr = oo), the base 
case, using the vorticity transport equation : 

^ = (w ■ V)« + i/V 2 « - w(V ■ u) + (1/^)(V ? x Vp(‘)) . (16) 

The four terms on the rhs of Eq. (16) describe, respectively, the rate of change of u? due to stretching 
and tilting of the vortex lines, viscous diffusion, velocity divergence due to density variation, and baroclinic 
torque. In variable density flows, this last term is nonzero only if the density gradient is not aligned with 
the pressure gradient. 

Since no mean shear is imposed on the present flow, both the vorticity and turbulence energy decay due to 
viscous action. In zero gravity, chemical reaction enhances the rate of vorticity decay in our flow for two 
reasons. First, the volumetric expansion, (V * u), due to the exothermic energy release is positive and thus 
acts as a sink in Eq. (16). Second, the magnitude of the baroclinic torque near the reaction zone is small 
since Vp and Vp are nearly aligned parallel to the vertical direction in our flow in which the constant- density 
surfaces and constant- pressure surfaces are nearly horizontal. The reduction of vorticity in a mixing layer 
due to non-premixed chemical reaction has been also reported by McMurtry et al. [6] and by Mahalingam 
et al . [5] in their study of forced, low speed, axisymmetric jet diffusion flame. 

The variance of the vorticity is the enstrophy which is completely determined from the turbulence energy 
spectrum, e.g. for incompressible flow the enstrophy is given by : 

w 2 (0/ 2 = D{t)/{ 2u)= [°° k 2 E(k,t)dk, (17) 

Jo 
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where D is the dissipation rate of the turbulence energy F, and k is the wave number. The transport 
equation for the enstrophy is 

D( ^ 2 / -- ) = (« • 5 • «) + i/V 2 (w 2 /2) - j/Vu, : (Vu,) T - w 2 (V • u) + (l//>> • (V#> x Vp«) , (18) 

Di 

where 5 is the rate of strain tensor. The terms on the rhs of Eq. (18) describe the rate of change of (w 2 /2) 
due to, respectively, vortex stretching or tilting, viscous diffusion, viscous dissipation, velocity divergence 
and baroclinic torque. It is seen that the velocity divergence and baroclinic torque have the same effect on 
the enstrophy as on the vorticity. Reducing the enstrophy means reducing the turbulence energy, E , and its 
dissipation rate, D. Reduction of D results in increasing the Kolmogorov scale, r;, and hence a reduction of 
scalar gradients, scalar mixing and thus the reaction rate [ 8 ]. 

Simulations for nonzero gravity with two values of Fr (17.6 and 7.04) were performed. 

Our results indicate that the baroclinic torque near the reaction zone increases with increasing the gravita- 
tional acceleration (i.e. reducing Fr). The evidence is clear both in the physical space results (contours of 
temperature, velocity and vorticity) and in the wave-number space (spectra of E(k) and D(k)). 

In the reaction zone, the buoyancy force causes hot products to rise in the 2 — direction and lower-temperature- 
fuel to move downwards toward the reaction zone. This is observed via the instantaneous contours of the w— 
velocity and temperature in any vertical ( xz ) plane. Thus alternating ‘hot’ and ‘cold’ fluid parcels cross the 
horizontal planes in the neighborhood of the reaction zone. A strong density gradient is thus established in 
the horizontal (x,y) directions in addition to that in the vertical direction. Accordingly Vp and Vp become 
misaligned thus increasing the magnitude of the baroclinic torque, the vorticity and the enstrophy. The 
instantaneous contours of u y show that it attains its maximum magnitude near the flame but it changes its 
sign across the reaction zone. The presence of the reaction zone between opposite-sign vorticities has been 
also observed by Mahalingam et al.[5] in their study of forced, low speed, buoyant axisymmetric jet diffusion 
flame. 

Now we discuss the effect of buoyancy on the spectra of turbulence energy. Since the reacting flow is inho- 
mogeneous in the z-direction, we can examine the spectra in horizontal planes. Figure 2 shows the spectra 
of the 2 -component of energy, £ 3 ( 6 ), at times t=2 and t=5 and x— component of energy, Fi(fc), at t=5, for 
Fr = 00 , 17.6, 7.04 on the plane 2 = 0.5 which is the initial interface between fuel and oxidant streams. It is 
seen that E$(k) decreases with time for the zero-gravity case, Fr = 00 , as explained earlier. In contrast, the 
energy and dissipation (not shown) increase with lowering Fr, thus reducing the Kolmogorov length scale 
and increasing the integral length scale, i.e. increasing the turbulence Reynolds number. The result is a 
reduction of the thickness of the reaction zone and an increase of its surface area due to wrinkling, hence an 
increase of the total consumption rate of the reactants. The increase of energy occurs first in E$(k) at high 
wave numbers (Fig. 2) where chemical reaction takes place, then due to thermal diffusion all scales become 
affected as time increases. The transfer of turbulence energy from E$(k) to both Ei(k) and F 2 W occurs 
via the pressure-strain interactions. The anistropy of turbulence is seen in Fig.2 by comparing E$(k) and 
Ei(k) at t=5. We also examined the balance of energy in the Fourier space: 

dE(k)/dt = T(k) + G(k) + P(k) - D(k) , (19) 

where the four terms on the rhs describe respectively the nonlinear transfer of energy in a given coordinate 
direction between the different scales, the transfer from potential to kinetic energy, the transfer of energy due 
to pressure gradients and the viscous dissipation. The time behavior of G(k) shows that initially high wave 
numbers (small scales) are the first to be affected by gravity, and in time the peak of G(k ) shifts gradually 
to lower wavenumbers (larger scales). Also G(k) is larger by at least an order of magnitude than the other 
terms on the rhs of Eq. (19) and increases monotonically as expected. 

4 Future Work 

The following phenomena will be studied in a stepwise manner during the next 3 years : 1) chemical kinetics 
for finite values of Da and examination of flame extinction in a range of values of Fr, Re\ ; 2) uniform mean 
shear ; 3 ) varying the direction of the gravity vector relative to the F/O interface plane (normal or parallel) 
; 4) dependence of fluid properties on temperature. 
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Figure 1: Schematic of flow configurations: (a) uniform mean velocity; (b) uniform mean velocity 
gradient in z— direction. The direction of the gravity vector may be either parallel or perpendicular 
to the initial interface between the fuel (F) and the oxidant (0), thus simulating a horizontal or 
vertical flame. 
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